Metallic clusters on a model surface: quantum versus geometric effects 



(N 

o 

(N 
+-> 

o 

O 



o3 



O 

o 



> 

^- 

O 

d 



X 



S. A. Blundell, 1 Soumyajyoti Haldar, 2,3 ' * and D. G. Kanhere 3 

l SPSMS, UMR-E CEA/UJF-Grenoble 1, INAC, Grenoble, F-38054, France. 
2 Centre for Modeling and Simulation, University of Pune, Ganeshkhind, Pune 411 007, India 
3 Department of Physics, University of Pune, Ganeshkhind, Pune 4.11 007, India. 

We determine the structure and melting behavior of supported metallic clusters using an ab initio 
density-functional-based treatment of intracluster interactions and an approximate treatment of 
the surface as an idealized smooth plane yielding an effective Lennard- Jones interaction with the 
ions of the cluster. We apply this model to determine the structure of sodium clusters containing 
from 4 to 22 atoms, treating the cluster-surface interaction strength as a variable parameter. For 
a strong cluster-surface interaction, the clusters form two-dimensional (2D) monolayer structures; 
comparisons with calculations of structure and dissociation energy performed with a classical Gupta 
interatomic potential show clearly the role of quantum shell effects in the metallic binding in this 
case, and evidence is presented that these shell effects correspond to those for a confined 2D electron 
gas. The thermodynamics and melting behavior of a supported Na2o cluster is considered in detail 
using the model for several cluster-surface interaction strengths. We find quantitative differences 
in the melting temperatures and caloric curve from density-functional and Gupta treatments of the 
valence electrons. A clear dimensional effect on the melting behavior is also demonstrated, with 2D 
structures showing melting temperatures above those of the bulk or (at very strong cluster-surface 
interactions) no clear meltinglike transition. 



PACS numbers: 65.80.-g, 61.46.Bc, 64.70. D-, 64.70. Nd 

I. INTRODUCTION 

Small particles and clusters can have markedly differ- 
ent physical and chemical properties from those of the 
bulk material because of the enhanced ratio of surface 
to volume and quantum confinement effects. [1] For ex- 
ample, the structure and atomic coordination, as well as 
electronic and magnetic properties, of small particles can 
all show new features. [2] In addition, the thermodynamic 
and melting behavior of clusters can have peculiar prop- 
erties. For instance, small clusters of Sn and Ga were 
found to undergo a meltinglike transition at tempera- 
tures higher than the melting point of the corresponding 
bulk material, [3, 4] contradicting the standard paradigm 
that a small particle should melt at a lower temperature 
than the bulk because of the effect of the surface. This 
behavior was explained in terms of a difference in the 
nature of the bonding between cluster and bulk, with 
small clusters of Sn and Ga both having a highly cova- 
lent character. [5, 6] Other experiments showed that small 
clusters of Na with from 50 to 300 atoms did undergo 
a mcltinglike transition at temperatures lower than the 
bulk melting point, but at temperatures that varied irreg- 
ularly with the size of the cluster, implying a competition 
between geometric and quantum effects. [7-11] 

Experimental[12-15] and theoretical[13, 14, 16, 17] 
work has also shown that the properties of clusters sup- 
ported on a surface can be modified by the interaction 
with the surface. However, there is a lack of work on the 
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thermodynamic properties and melting behavior of sup- 
ported clusters using first-principles, density-functional- 
based electronic structure methods. 

A recent motivation for studying theoretically the 
structural and thermal properties of supported metal- 
lic clusters is provided by the cluster-catalyzed growth 
process of carbon nanotubes (CNTs).[18, 19] In this pro- 
cess, the CNT grows from a small cluster, typically of a 
transition metal (or its oxide), supported on a substrate. 
It is highly desirable to be able to produce large quanti- 
ties of high-quality CNTs with controlled properties, [18] 
but many aspects of the growth process are poorly un- 
derstood. For instance, the efficiency of the growth pro- 
cess can depend markedly on the particular combination 
of substrate and cluster material chosen, as well as on 
other experimental parameters such as the temperature 
and pressure. Much recent theoretical work has been 
devoted to the simulation of the cluster-catalyzed CNT 
growth process. [20-25] 

In this context, Ding et al. [22] and Shibuta and 
Maruyama[23] considered a simplified model in which the 
surface was replaced by an idealized smooth plane inter- 
acting with the cluster via a potential of Lennard- Jones 
type. Treating the cluster-surface interaction (Lennard- 
Jones well depth) as a variable parameter, they showed 
that the melting temperature of clusters with several 
hundred atoms increased monotonically as the strength 
of the cluster-surface interaction increased. Using a sim- 
ilar model, Sarkar and Blundell[26] considered smaller 
clusters of size 55 atoms and showed that there were 
detailed changes in the structure of the cluster as the 
cluster-surface interaction varied, accompanied by steps 
in the melting temperatures. Also, Jiang et oL,[25] in 
a study of the thermodynamics of supported Fe-C clus- 
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ters, have employed a smooth plane interacting with the 
cluster via an effective Morse potential, the parameters 
of which were fitted to ab initio calculations. 

Now, the above-mentioned calculations[22, 23, 25, 26] 
all use classical molecular dynamics (MD)[27] with para- 
metric interatomic potentials to simulate the thermo- 
dynamic properties of the supported clusters. In this 
work, we retain the simplified model of the surface, 
but employ instead an ab initio treatment of the in- 
tracluster interactions, within density functional theory 
(DFT),[28] directly in both the thermodynamic simula- 
tions and the global optimization procedure used to de- 
termine the lowest-energy cluster structure at K. The 
cluster-surface interaction is treated as an effective clas- 
sical external potential acting on the ions of the cluster. 
While DFT-based thermodynamic simulations are orders 
of magnitude more expensive than classical MD simula- 
tions, they are nevertheless feasible and have been used 
with success in the past, for example, to explain the pe- 
culiar melting behavior of free (unsupported) Sn and Ga 
clusters. [5, 6] Also, first principles determinations of the 
melting temperature of small Na clusters were made in 
Refs. 29 and 30. To illustrate our general approach, in 
this work we revisit Na clusters, with the aim of show- 
ing that the model is capable of accounting for subtle 
quantum effects in the metallic bonding. 

The plan of the paper is as follows. In the next sec- 
tion we describe in some detail our model and calculation 
procedures. Then, in Sec. Ill, we consider the particular 
case of a cluster-surface interaction that is sufficiently 
strong that the clusters collapse into a monolayer, or 
two-dimensional (2D), structure. This gives a particu- 
larly illuminating example of the differences between a 
classical and a quantum treatment of the valence elec- 
tron gas. For this case, we determine the lowest-energy 
structures for sizes 4 < N < 22. In order to bring out the 
quantum effects in the metallic bonding, we also perform 
calculations for comparison with a classical interatomic 
potential. In addition to revealing a competition between 
geometric and quantum effects, in this section we also 
explore the properties of '2D metallic clusters' and con- 
sider evidence for 2D (rather than 3D) quantum effects. 
Next, in Sec. IV, we consider in detail the thermody- 
namics and melting behavior of a supported Na2o cluster 
within the model for three values of the cluster-surface 
interaction strength, once again bringing out the similar- 
ities and differences with the same simulations performed 
with classical interatomic potentials. We also explore the 
dimensional effect on the melting behavior when the clus- 
ters are predominantly 2D. The conclusions are given in 
Sec. V. 



II. MODEL AND METHODOLOGY 

A simple system displaying the main physical proper- 
ties of metallic clusters is provided by Na, [1] and we shall 
consider Na clusters throughout. Our main approach 



involves a DFT-based Kohn-Sham (KS) description [28] 
of the metallic cluster. We employ either the 
local-density approximation (LDA) or the generalized- 
gradient approximation (GGA) with Vanderbilt's ultra- 
soft pseudopotential,[31] as implemented in the vasp 
package. [32] These approaches will be denoted KS-LDA 
and KS-GGA, respectively, in the following. In parts 
of the work, we also use a simplified (and computation- 
ally faster) real-space KS method in the LDA incorpo- 
rating a soft, phenomenological, local pseudopotential 
of the form described by Blaise et al. [33] This approx- 
imate KS method (henceforth referred to as 'KS-soft') 
has been used with success to describe the fragmenta- 
tion of charged Na clusters [34] and the melting of free 
(unsupported) Na clusters. [35] The KS formalism yields 
an expression [2 8] for the internal energy E c \ us of the free 
cluster (that is, not yet taking into account explicitly the 
interaction with the surface) . 

For comparison, we also describe the metallic bond- 
ing within a Na cluster by a classical many-body Gupta 
potential derived within the second-moment approxima- 
tion (SMA) . [36-38] The internal energy of a free (unsup- 
ported) cluster in this approach is given by 
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where Rij is the distance between ions i and j. The first 
term in Eq. (1) is a repulsive potential of Born-Mayer 
form, and the second term is a cohesive energy. We take 
the parameters for Na from the work of Li et al.:[39] 
£ = 21.398 mRy, A = 1.1727 mRy, p = 10.13, q = 1.30, 
and ro = 6.99 ao, where ao is the Bohr radius. These 
authors fitted the parameters to a database of total en- 
ergies as a function of lattice constant (for both fee and 
bec lattices) obtained from calculations in the LDA for 
bulk Na. The potential was then checked by using it 
to predict bulk properties such as the equilibrium lattice 
constant and bulk modulus. For instance, the bulk melt- 
ing temperature T m predicted by the model was found to 
be 333 K, compared to the experimental value of 371 K. 

Following Ding et al. [22] (and the general approach of 
Shibuta and Maruyama[23]), the surface is modeled as 
an idealized smooth plane that interacts with the cluster 
via a Lennard- Jones 9/3 potential, yielding an interaction 
energy 
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where Zi is the coordinate of ion i perpendicular to the 
surface. This model for the surface is used with both the 
SMA (1) and the KS treatments of the intracluster inter- 
actions used to determine E c i us . In this way, the cluster 
is constrained to lie in the vicinity of the minimum of 
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the Lennard-Jones well, which is located roughly a dis- 
tance a above the Z = plane. We choose the param- 
eter a = 0.3 nm. As was done in Refs. 22 and 23, the 
Lennard-Jones well depth s may be treated as a variable 
parameter describing the overall strength of the cluster- 
surface interaction. The total energy of the cluster plus 
surface is 

Etot — E c \ us + E[ nt . (3) 

Using these models of a supported metallic cluster, we 
carry out a global optimization procedure to search for 
the structure that minimizes E tot , and perform statisti- 
cal simulations to extract thermodynamic averages and 
study the melting behavior of the cluster. With both 
the SMA potential (1) and the KS-soft method, [34, 35] 
we determine the optimum structure at K by means of 
a basin- hopping algorithm. [40] This algorithm involves a 
Monte-Carlo (MC) process, combined with optimization 
to the nearest local minimum at each MC step, to ex- 
plore the minima on the potential-energy surface. With 
the SMA potential, we use 10 5 -10 6 hops (local minimiza- 
tion steps). The basin- hopping procedure generates a 
sampling of excited isomers as well. We store the first 
(i.e., lowest-energy) six excited isomers found with the 
SMA potential in a library of structures, together with 
four more higher-energy isomers chosen randomly from 
the complete set of local minima generated. 

Now, the KS model is orders of magnitude more expen- 
sive computationally than the SMA model, and to find 
the optimum cluster structure within the KS approach 
we proceed as follows. Using the library of structures 
generated with the SMA potential as seeds (initial struc- 
tures) for the MC process, we perform 50-100 basin hops 
within the KS-soft model for each seed structure. Thus, 
we perform in total 500-1000 basin hops for each cluster 
size N, the smaller number applying to the larger clus- 
ters that we consider having up to N = 22 atoms. For 
the smaller cluster sizes N < 10, we find that on dou- 
bling the total number of basin hops, we do not observe 
a change in the lowest-energy structure found. Note that 
a basin-hopping process based directly on the KS model, 
as we perform here, often yields structures that were not 
present at all among the isomers of the SMA model (see 
Sec. IV for examples), and therefore we do not simply 
relax large numbers of isomers found within the SMA 
model. To speed up the basin-hopping process, we use 
a relatively small simulation box for the KS solution of 
side 45 ao- For the larger sizes N > 8 considered, there 
are usually many isomers with closely-spaced total ener- 
gies (see Sec. IV), and in some cases the ordering of these 
isomers can be modified by confinement effects for a sim- 
ulation box this small. Therefore, we perform a second 
step in which the ten lowest-energy structures found from 
the basin-hopping step are relaxed to a precise energy tol- 
erance (6 E to t ~ 10 -8 Ha) using a larger simulation box 
of side 90 ao (but the same real-space grid spacing). At 
this stage, in some instances (as described in the text) 
we also relax the structures found within the KS-LDA or 



KS-GGA treatments. [32] 

Finally, we perform thermodynamic simulations for 
Na 2 o following the general procedures described in our 
earlier work, [35] using both the SMA and KS-LDA mod- 
els of the cluster. With the SMA model, we perform 
of order 30 constant-total-energy (microcanonical) MD 
simulations distributed over a range of kinetic tempera- 
tures from 30 K to around 750 K, with each simulation 
of order 1 ns (so that the total simulation time is around 
30 ns per cluster). For kinetic temperatures higher than 
about 700 K, the Na clusters tend to evaporate on this 
time scale. We start at low kinetic temperature using the 
optimum structure found above, and increase the kinetic 
temperature gradually from run to run, using the coor- 
dinates and rescaled velocities at the end of one run to 
provide the initial condition for the next run. Thus, in ef- 
fect we slowly heat the cluster. The kinetic temperatures 
chosen for the MD simulations are more closely spaced in 
the range 200 K to 450 K, where the cluster meltinglikc 
process tends to occur. After this, we perform a multiple 
histogram fit to the overlapping histograms of potential 
energy from the various simulations to extract the classi- 
cal ionic density of states il(E). With fl(E) in hand, one 
can now evaluate thermodynamic averages such as the 
ionic specific heat in a variety of ensembles, including 
the canonical ensemble. [35] 

In the case of the KS-LDA cluster model, we also 
proceed via a slow-heating algorithm, this time per- 
forming a sequence of isokinetic Born-Oppenhcimer MD 
simulations [41] at gradually increasing kinetic energies. 
Because the KS-LDA approach is much more expensive 
than SMA, we use a total simulation time per cluster that 
is somewhat smaller than for the SMA model, of about 
1 ns to 3 ns. Finally, a canonical multiple histogram fit 
is used to extract fi(i£).[35] 



III. MONOLAYER STRUCTURES FOR 

N = 4-22 

As we shall see in the next section, when the cluster- 
surface interaction e [see Eq. (2)] is increased, the clus- 
ter tends to become progressively flatter. Useful insight 
into the role of quantum versus geometric effects can be 
gained by considering a cluster-surface interaction suffi- 
ciently strong that the lowest-energy structures are all 
monolayer. For this purpose we have chosen e = 0.5 eV. 
We have found optimum structures for both the SMA 
model of intracluster interactions, Eq. (1), and the KS- 
soft model[34, 35] using the basin-hopping algorithm, as 
described in Sec. II. The final structures are shown in 
Fig. 1. We also give the dissociation energies A£?diss (the 
energy required to remove a single neutral atom from the 
cluster) for these structures in Fig. 2. Noting that the 
atoms of a monolayer structure all lie at the minimum — e 
of the Lennard-Jones well, Eq. (2), we have here defined 

AE diss (N) = E tot (N - 1) - Et t(N) -e, (4) 
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correcting for the trivial contribution arising from the 
cluster-surface interaction. Thus I\E^ SS (N) represents 
just the 'internal' contribution to the dissociation energy 
due to the intracluster interactions. 

We see that the SMA structures in Fig. 1 form a geo- 
metrical packing in a triangular lattice, favoring in par- 
ticular hexagonal cluster structures. Thus, for N = 7 
and N = 19 atoms the structures are two- and three-shell 



hexagonal arrangements, respectively. These are exam- 
ples of compact, 'closed-geometric-shell' structures in 2D. 
Similarly, the structure for N = 10 is a combination of 
two hexagonal structures. The structures either side, for 
example those at N = 18 and N — 20, are formed by 
adding or subtracting one atom to or from the outside 
these closed-shell geometric structures. In fact, a strik- 
ing property of the SMA structures in Fig. 1 is that for 
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FIG. 2. (Color online) Dissociation energies A_Bdi ss of the 
supported monolayer structures of Fig. 1, corrected for the 
cluster-surface interaction strength e = 0.5 eV. 



all 4 < N < 21 the structure for any size N is obtained 
simply by adding one atom to the outside of the struc- 
ture for N — 1. This simple growth sequence is broken 
for the first time at N = 22. 

The dissociation energy of the SMA structures (see 
Fig. 2) mirrors this growth sequence. Thus, for the 
smaller sizes there are peaks at N = 7 and 10, re- 
flecting the stability of these compact geometric struc- 
tures. There is then a sawtooth pattern in AE^^N) for 
N = 10-18 and N = 19-22. We can understand this saw- 
tooth variation in a very simple way by counting the num- 
ber of nearest-neighbour bonds in Fig. 1 broken in each 
dissociation. For N — 10-18, the atom to be removed 
from an even-iV structure is connected by three bonds 
to the remainder of the cluster, while that for an odd-iV 
structure is connected by only two. In consequence, the 
even- TV structure has a higher dissociation energy. The 
sawtooth is interrupted at N — 18 and 19; because of the 
geometrical properties of the growth sequence, the atom 
to be dissociated has three bonds for both of these sizes. 
The sawtooth pattern then resumes after N = 19, but 
this time it is the odd sizes N that have three bonds and 
the greater dissociation energy. Note also that it tends 
to be an atom on a corner that is removed in a dissoci- 
ation. This happens because it is these atoms that have 
two or three bonds to the rest of the cluster; an atom in 
the center of a complete side has four bonds. 

Turning to the KS structures in Fig. 1, we note that 
there are some similarities and some differences with the 
SMA structures. Up to size N = 13, the KS structures 
also follow a simple growth pattern in which atoms are 
added successively to the outside of each structure, al- 
though the structures for N = 6 and 9 differ from their 
SMA counterparts. For larger sizes iV > 13, this simple 
growth sequence is broken. Some larger structures coin- 
cide with their SMA counterpart (N = 19 and 21), but 
the others generally prefer a more elongated form. Also, 
the KS dissociation energies in Fig. 2 form a consistent 
sawtooth pattern for all 4 < N < 22. 



These differences between the KS and SMA approaches 
are indicative of the role of quantum effects in the system 
of valence electrons that forms the metallic bonding in 
the cluster. The SMA potential favors a purely geometri- 
cal packing, and the dissociation energies are sensitive to 
the number of nearest-neighbor bonds. In metallic clus- 
ters, however, the valence electrons display a fermionic 
'shell' structure — a finite-size quantum effect — which is 
particularly pronounced at the smaller sizes considered 
here.[l, 2] A cluster with a closed fermionic shell is par- 
ticularly stable, in analogy to the noble gases in the peri- 
odic table. In the case of free (unsupported) Na clusters, 
the closed shells coincide with those of the 3D simple 
harmonic oscillator (SHO), namely, for N = 2, 8, and 
20 electrons. [1, 2] Since in a Na cluster each atom con- 
tributes one valence electron, these are also the ('magic') 
numbers of atoms for which there is pronounced stability. 
For the monolayer structures discussed here, one might 
expect that it is the closed shells of the 2D SHO that are 
appropriate: N — 2, 6, 12, and 20. These magic numbers 
occur in quasi-2D semiconductor quantum dots. [42] 

For a closed-shell system, the valence-electron gas is 
generally circular (in 2D) or spherical (in 3D), but be- 
tween closed shells, the valence-electron gas tends to 
minimize its energy by spontaneously deforming. This 
phenomenon is well-known in nuclear physics, [43] and 
has been shown to apply also to small, free 3D metal- 
lic clusters. [44] For small cluster sizes, the deformation 
of the valence-electron gas can in turn drive the distri- 
bution of Na + ions away from circular symmetry, pro- 
ducing not just a distortion, but a new structure. We 
see evidence for this phenomenon in Fig. 1 in the ten- 
dency of the KS model to yield elongated ground-state 
structures (relative to the SMA structures) for N > 12, 
as is especially noticeable in the size range N = 14-18. 
(The case N = 20, discussed in more detail in the next 
section, is an exception, because N = 20 is expected to 
have closed fermionic shells, yet the optimum KS struc- 
ture is still somewhat elongated.) We also see from Fig. 1 
that the compact geometric structure can still have the 
lower energy in particularly favorable cases, for example, 
the three-shell hexagonal structure at N = 19. Thus, 
there is evidence of competition between geometric and 
quantum effects in the KS structures. 

The role of quantum effects is particularly striking in 
the dissociation energy, Fig. 2, where one finds a consis- 
tent sawtooth pattern in which the odd sizes always have 
the smaller dissociation energy. This effect, observed also 
in free metallic clusters, [1, 2] is a result of the pairing of 
spins: to a first approximation, the spins for N even form 
pairs (spin up and down) in each single-particle energy 
level, yielding a more stable structure than for N odd, 
which has an unpaired final spin. Note that while the 
SMA model also yielded a sawtooth pattern for some 
sizes, this had a geometric origin and it could be either 
N even or N odd that gave the smaller dissociation en- 
ergy. Also, the amplitude of the odd-even oscillations in 
Fig. 2 tends to be greater for the KS model. 



There may also be indications of 2D fermionic shell 
closures in the KS dissociation energies in Fig. 2. It is 
striking that the amplitude of the odd-even oscillations 
decreases abruptly for N > 12, which could be because 
a new fermionic shell has been begun following a closed 
shell at N = 12. The evidence for a 2D fermionic shell 
at N = 6 is less clear, although the dissociation energy 
for the system with one additional atom, N = 7, is par- 
ticularly low. The evidence for a shell closure at N = 20 
seems unclear, however, either in the structures or in the 
dissociation energy. 



IV. THERMODYNAMICS AND MELTING FOR 
N = 20 

In this section, we discuss in detail the structure and 
thermodynamics of a supported Na2o cluster within our 
model. Figure 3 shows the lowest-energy structures 
of Na2o as a function of the cluster-surface interaction 
strength e [see Eq. (2)]. The intracluster interactions 
in this figure are treated within the KS-LDA approach. 
After determining candidate lowest-energy structures at 
values of e from to 0.5 eV in steps of 0.05 eV, we made 
a final series of high-precision relaxations in order to plot 
the total energy (3) of each structure found as a func- 
tion of e (see the upper panel of Fig. 3). The crossovers 
of these energy curves enable us to identify a range of 
e values for which a particular structure is the lowest- 
energy structure (that we have found). In this way, we 
find a series of structural transitions at particular values 
of e, in which the cluster becomes progressively flatter 
as e increases. This behavior is similar to that found for 
small clusters of Fe, Co, and Ni in Ref. 26 using classical 
methods with parametric interatomic potentials. 

The ground-state geometry of free Na2o (e = 0), 
Fig. 3(a), is quite spherical in shape. For small values 
of the cluster-surface interaction strength e « 0.02 eV, 
as the cluster-surface interaction starts to become com- 
parable to the internal interactions within the cluster, 
the lowest-energy structure deforms from the spherical 
shape to a slightly elongated, though still quite com- 
pact, form with the largest face aligned parallel to the 
surface [Fig. 3(b)]. As e increases further, the lowest- 
energy structure eventually flattens, first into a structure 
with three ionic layers [Fig. 3(c)], and then into a two- 
layer structure [Fig. 3(d)]. Finally, the structure becomes 
monolayer for e > 0.36 eV [Fig. 3(e)]. The flatter struc- 
tures lower their energy by increasing the contact area 
between the cluster and the surface for larger e. 

Let us consider the isomers in more detail in one case, 
that of the monolayer structure, Fig. 3(e). Figure 4 shows 
the ground and first three excited isomers in the classical 
SMA potential and several KS models. The SMA isomers 
have been determined after 10 6 basin hops. For the KS 
isomers, we first performed of order 500 basin hops within 
the KS-soft model (as discussed in Sec. II), retaining the 
lowest seven structures found. A final high-precision re- 
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FIG. 3. (Color online) Total energy of a Na2o cluster on a 
surface (upper panel) as a function of cluster-surface interac- 
tion strength e [see Eq. (2)] for the structures (a)-(e) shown 
in the lower panel. Structure (a) (short-dashed line) is the 
lowest-energy structure for < e < 0.02 eV; structure (b) 
(full line) for 0.02 eV < e < 0.18 eV; structure (c) (dotted 
line) for 0.18 eV < £ < 0.21 eV; structure (d) (full line) for 
0.21 eV < e < 0.36 eV; and structure (e) (dashed line) for 
e > 0.36 eV. 

laxation of these structures was made in the KS-LDA and 
KS-GGA approaches [with a projected-augmented plane- 
wave (PAW) treatment of ionic cores, as implemented in 
vasp[32]]. 

The three lowest-energy structures within the SMA 
[Figs. 4(a)-(c)] are nearly degenerate to within about 
1 meV, reflecting the fact that these structures may be 
regarded as surface rearrangements of each other. Thus, 
it is possible to obtain the ground-state SMA structure 
for N = 21 in Fig. 1 by adding one atom to the outer 
shell of each of these three structures. 

The situation is different with a KS treatment of in- 
tracluster interactions. The KS ground-state structure 
[Fig. 4(e)] corresponds to an excited isomer within SMA 
(lying at 20 meV), and is separated from the first excited 
isomer by about 20-35 meV. (While there is some varia- 
tion in the excitation energy of isomers according to the 
KS model used, the ground-state structure is the same 
for each model.) The first excited isomer in the SMA 
[Fig. 4(b)] appears to be disfavored in the KS approach; 
it is not found by the basin-hopping procedure, and when 
we tried to relax it within a KS approach it 'slipped' to 
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FIG. 4. (Color online) Ground-state (left column) and first 
three excited (right columns) monolayer isomers of Na2o 
found by the basin-hopping procedure for a cluster-surface 
interaction strength e = 0.5 eV. Notation: SMA, second- 
moment approximation; KS, Kohn-Sham; PAW, projected- 
augmented plane wave; GGA, generalized gradient approxi- 
mation; LDA, local density approximation; KS-soft, the sim- 
plified Kohn-Sham model in the LDA with a soft, phenomeno- 
logical pseudopotential (see text). Excitation energies AE are 
shown relative to the ground-state structure. 



the geometry shown as the first excited isomer for KS 
[Fig. 4(f)]. Note that the first three excited isomers in 
the KS model [Figs. 4(f)-(h)] may be regarded as sur- 
face rearrangements of each other and of the lowest three 
SMA isomers [Figs. 4(a)-(c)]. However, the separation of 
these three isomers of KS is much greater than the 1 meV 
separation found for the SMA structures. These obser- 
vations highlight the fact that the energy of the metallic 
cluster in a quantum approach is sensitive, in a compli- 
cated way, to the deformation and symmetry of the wave 
function, and not just to the geometric packing of the 
ions or the number of nearest neighbors. 

Note that excited structures such as (g) in Fig. 4, 
or ground-state structures such as the KS structure for 
N = 16 in Fig. 1, are not found among the structures 
given by the SMA basin-hopping process (even after 10 5 - 
10 6 hops), but become common in a KS basin- hopping 
process. This illustrates the importance of basing the 
basin-hopping process directly on the KS model in order 
to achieve a fully unbiased search for the ground-state 
structure, despite the high computational cost of such an 
approach. 

Finally, we turn to the thermodynamic and melting 
behavior of the supported Na2o cluster. Figure 5 shows 
the canonical ionic specific-heat capacity of the cluster 
for e = 0.1 eV, extracted by an ab initio MD simulation 
and multiple-histogram fit using a KS-LDA description 
of the intracluster interactions throughout (as described 
in Sec. II). For comparison, we also give the specific- heat 
curve for a free Na2o cluster (e = 0) calculated by similar 
methods, taken from Ref. 30. Inspection of the ionic tra- 
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FIG. 5. (Color online) Canonical ionic specific-heat capacities 
of Na2o, for surface-cluster interaction strength [see Eq. (2)] 
e = 0.1 eV and for the free cluster (e = 0, taken from Ref. 30), 
for KS descriptions of the intracluster interactions. The spe- 
cific heat is expressed as a multiple of its zero-temperature 
(classical) limit, Co- 



jectories shows a meltinglike transition to occur for both e 
values, with the cluster passing from a solidlike behavior 
at low temperature (vibration of ions around fixed points 
combined with overall rotation) to a liquidlikc behavior 
at high temperatures (diffusion throughout the entire vol- 
ume of the cluster). The meltinglike transition is broad, 
with a width of around 50-100 K and a 'melting temper- 
ature' (conventionally corresponding to the maximum on 
the specific-heat curve) around 220-240 K. There is a sig- 
nificant change in the specific-heat curve between the free 
and supported cluster: for e = 0.1 eV the curve is broader 
and the meltinglike transition correspondingly less well- 
defined, and the melting temperature shifts to a slightly 
higher value. This marked change in the specific-heat 
curve is associated with a change in the lowest-energy 
structure at K (see Fig. 3). 

Canonical ionic specific-heat curves are shown in 
Figs. 6 and 7 for e = 0.38 eV and e = 0.5 eV, respectively. 
For both these values of e, the lowest-energy structure at 
K is a monolayer structure (see Fig. 3), but the choice 
e = 0.38 eV for Fig. 6 is very close to the critical value 
where the lowest-energy structure becomes a two-layer 
structure. This leads to an interesting difference in the 
melting behavior. 

For e — 0.38 eV, inspection of the ionic trajectories 
shows that, when the cluster is liquidlike, ions frequently 
hop to the second layer and then back to the first, usu- 
ally reentering the first layer at a different position from 
where they left. This process occurs readily because there 
are energetically close isomers with a two-layer structure. 
We find a clear peak in the specific-heat capacity around 
480 K, corresponding to the temperature at which the 
ions first have sufficient energy (on average) to hop to the 
second layer (that is, to overcome the barrier required to 
move to an energetically close two-layer structure). For 
higher temperatures, the cluster is found to be in a liq- 
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FIG. 6. (Color online) Canonical ionic specific heat capacity 
of Na2o, with surface-cluster interaction strength [see Eq. (2)] 
e = 0.38 eV, for SMA and KS descriptions of the intracluster 
interactions. The specific heat is expressed as a multiple of 
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FIG. 7. (Color online) Canonical ionic specific heat capacity 
of Na2o, with surface-cluster interaction strength e = 0.5 eV. 
See Fig. 6 for other notation. 



orate on the time scale of several hundred picoseconds 
used for the sampling runs. 

We thus demonstrate clear dimensional effects in the 
melting of the clusters. When the ions are constrained 
to move mostly in 2D, as they are at e = 0.5 eV, the 
possibilities for diffusion are reduced compared to the 
3D case (free Na clusters), and the onset of a liquidlikc 
behavior is gradual, with no obvious peak in the specific- 
heat curve. At intermediate values of e « 0.4 eV, the 
motion is still predominantly 2D, but there are energet- 
ically close two-layer isomers; when the barrier to these 
configurations can be overcome, the sudden availability 
of additional phase space can give a peak in the specific- 
heat curve and a meltinglike transition. Even in this case, 
however, we note that the observed melting temperature 
of T rn 480 K is significantly above the melting point 
of bulk Na, 371 K, contrary to the case of small free Na 
clusters, [7, 8] which melt at temperatures below that of 
the bulk (and contrary to the general paradigm for the 
melting of a small particle [1]). 

For comparison, we also give in Figs. 6 and 7 the spe- 
cific heat calculated for SMA interatomic interactions. 
Inspection of the ionic trajectories shows that the mech- 
anisms of melting are qualitatively similar to the mecha- 
nisms in the KS model, but there are quantitative differ- 
ences in the specific-heat curves. Thus, for e = 0.38 eV 
the meltinglike transition occurs around 260 K, a lower 
temperature than observed with the KS model. Similar 
differences in melting temperatures were observed for free 
(unsupported) Na clusters in Refs. 29 and 30. In these 
references it was found that ab initio DFT-based melting 
temperatures were in significantly better agreement with 
experiment than those calculated with the SMA poten- 
tial. 



V. CONCLUSIONS 



uidlike state, and the mechanism of hopping to the sec- 
ond layer and back contributes significantly to the overall 
diffusion of the ions within the cluster. 

On the other hand, for e = 0.5 eV, one finds that hops 
to the second layer are highly suppressed. At 300 K, for 
instance, the system spends less than 1% of its time in 
a two-layer configuration, and the ionic dynamics is thus 
constrained mostly to a plane (2D). At the lower temper- 
atures T < 300 K, parts of the ground-state triangular 
lattice become distorted owing to the random movement 
of the ions, yielding distorted triangles or rectangles, for 
example. The ionic motion becomes gradually more dif- 
fusive as the temperature increases. For temperatures 
above 400 K, one finds that the ions also begin to move 
increasingly in a vertical direction, toward the second 
layer. The increasing phase space thus available to the 
ions means that the specific heat continues to rise beyond 
T = 400 K, and there is no clear peak in the specific heat. 
Around T > 600 K, the ions of the cluster tend to evap- 



We have presented a model for a supported metallic 
cluster that employs a first principles DFT-based descrip- 
tion of the cluster, but approximates the surface as an 
idealized smooth plane interacting with the cluster via a 
parametric potential (here taken to be of Lennard- Jones 
type). Thermodynamic simulations, or global optimiza- 
tion methods such as basin hopping, using this model 
are comparable in computational cost to similar meth- 
ods with the free clusters. The results of this model were 
compared in detail with one in which the intracluster in- 
teractions were described instead by a parametric inter- 
atomic potential of SMA type, showing that important 
quantum effects were absent in the SMA-type simulations 
for small metallic clusters. The classical SMA model fa- 
vors a purely geometrical packing and is sensitive to the 
number of nearest neighbors, while a quantum treatment 
of the valence electrons (as provided by a KS model) 
yielded significant differences both in the lowest-energy 
structures and in energetic properties such as the dissoci- 
ation energy of the supported clusters. These differences 



9 



could be ascribed in part to fcrmionic shell closures, and 
we found some evidence that these shell closures could be 
2D rather than 3D for the monolayer cluster structures 
that occur for large cluster-substrate interactions. 

We showed that the Kohn-Sham and Gupta models 
could yield significant quantitative differences in melting 
temperatures and caloric curves. We also demonstrated 
clear dimensional effects in the melting of supported clus- 
ters by considering the case of a high cluster-surface in- 
teraction strength yielding a 2D monolayer structure. 
For intermediate values of the cluster-surface interaction 
strength e w 0.4 eV, the Na2o cluster gave a clear melt- 
inglike transition, with a peak in the specific-heat curve, 
but at an above-bulk melting temperature. In this case, 
there were energetically close two-layer cluster struc- 
tures. At higher cluster-surface interaction strengths 
e > 0.5 cV, where the ionic motion is constrained to 
be highly 2D, the onset of liquidlike behavior is grad- 
ual and there is no clearly definable 'melting tempera- 
ture.' Even for small values of the cluster-surface in- 
teraction strength, e ss 0.1 eV, the melting temperatures 
and caloric curves could be markedly different from those 
of the free cluster. 

An obvious generalization of the present model would 
be to consider an atomistic surface, and to extend the 
DFT treatment to the atoms of the surface. Such an ap- 



proach for thermodynamic properties would be computa- 
tionally very expensive, however, and the present model 
is a compromise solution, useful for a class of systems. It 
is likely to be realistic whenever the metallic bonding is 
not significantly disrupted by interaction with the sub- 
strate, which may be the case for an insulating substrate, 
for example. Our approach might also be useful in sit- 
uations where the main interest is in the properties of 
the cluster away from the surface, such as the interaction 
between the atoms of a carbon nanotube and the clus- 
ter in a description of the nanotube growth process, [20] 
where the nanotube grows from the top of the supported 
cluster. 
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